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ABSTRACT 

The powerful high-energy phenomena typically encountered in astrophysics in- 
variably involve physical engines, like neutron stars and black hole accretion disks, 
characterized by a combination of highly magnetized plasmas, strong gravitational 
fields, and relativistic motions. In recent years numerical schemes for General Rela- 
tivistic MHD (GRMHD) have been developed to model the multidimensional dynamics 
of such systems, including the possibility of an evolving spacetime. Such schemes have 
been also extended beyond the ideal limit including the effects of resistivity, in an at- 
tempt to model dissipative physical processes acting on small scales (sub-grid effects) 
over the global dynamics. Along the same lines, the magnetic field could be ampli- 
fied by the presence of turbulent dynamo processes, as often invoked to explain the 
high values of magnetization required in accretion disks and neutron stars. Here we 
present, for the first time, a further extension to include the possibility of a mean-field 
dynamo action within the framework of numerical 3 + 1 (resistive) GRMHD. A fully 
covariant dynamo closure is proposed, in analogy with the classical theory, assuming 
a simple a-effect in the comoving frame. Its implementation into a finite-difference 
scheme for GRMHD in dynamical spacetimes [the X-ECHO code: Bucciantini & Del 
Zanna (2011)] is described, and a set of numerical test is presented and compared with 
analytical solutions wherever possible. 

Key words: methods: numerical - (magnetohydrodynamics) MHD - magnetic fields 
- relativity - gravitation. 



1 INTRODUCTION 

A strong magnetic field plays a crucial role in many high- 
energy astrophysical systems. It is believed to be the key ele- 
ment in the context of Gamma Ray Bursts (GRBs) (Duncan 
& Thompson 1992; Thompson 1994; Meszaros & Rees 1997; 
Lee, Wijers & Brown 2000; Lyutikov & Blackman 2001; Vla- 
hakis & Konigl 2001; van Putten & Levinson 2003; Lyu- 
tikov 2006b; Komissarov & Barkov 2007; Metzger, Thomp- 
son & Quataert 2007; Uzdensky & MacFadyen 2007a,b; 
Barkov & Komissarov 2008; Bucciantini et al. 2008, 2009; 
Lyons et al. 2010; Metzger et al. 2011; Rezzolla et al. 2011), 
AGN-jets (Blandford & Znajek 1977; Blandford & Payne 
1982; Belvedere & Molteni 1984; Khanna & Camenzind 
1992; Konigl & Kartje 1994; Balbus & Hawley 1998; Tomi- 
matsu 2000; Fendt & Memola 2001; Sauty, Tsinganos & 
Trussoni 2002; van Putten & Levinson 2003; Duttan & Bier- 
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mann 2007; Hawley 2008; Penna et al. 2010; Tchekhovskoy, 
Narayan & McKinney 2011), magnetars (Duncan & Thomp- 
son 1992; Thompson & Duncan 1996; Murakami 1999; Lyu- 
tikov 2003, 2006a; Braithwaite & Spruit 2006; Woosley 2010; 
Kasen & Bildsten 2010; Yu 2011), and its dissipation and re- 
connection to be at the origin of many typical high energy 
phenomena (Uzdensky 2011). A large-scale ordered mag- 
netic field is in fact crucial to power many high energy sys- 
tems. A strong magnetic field close to the central black hole 
(BH) in accretion disks is invoked to explain the launching 
of relativistic jets. The pulsar and magnetar paradigms re- 
quire a strong dipolar magnetic field in neutron stars (NSs) . 
Magnetic field's stresses provide an efficient way to convert 
rotational or accretion energy into bulk flow, and to power 
relativistic winds. 

However, the origin of this strong and ordered magnetic 
field remains poorly understood. In particular, the environ- 
ment in which such a magnetic field is supposed to arise 
is often characterized by turbulence (Zhang, MacFadyen & 
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Wang 2009; Mizuno et al. 2011) and instabilities, ranging 
from MRI in accretion disks (Balbus & Hawley 1998) to 
kink and Tayler and convection in NSs (Miralles, Pons & 
Urpin 2000, 2002; Braithwaite & Spruit 2006; Ruediger et al. 
2008). While turbulent small-scale motions can easily am- 
plify a seed field up to equipartition with the turbulent ki- 
netic energy, one expects such a field to be highly tangled. If 
a large-scale ordered mean field can arise due to small scale 
velocity fluctuations is highly debated, and understanding 
what its configuration and its geometry are, and under what 
conditions it is realized, is of fundamental importance. It is 
not unreasonable that, during the formation of compact ob- 
jects like BHs or NSs, any frozen-in large-scale field might 
be amplified due to advection. However, as an explanation 
of the required magnetic field, this simply shifts the problem 
from the object to the progenitor. Moreover, in BHs, NSs, 
and GRBs, a large amount of angular momentum is required 
in the engine, but any strong and pre-existing magnetic 
field would rapidly slow down rotation (Duez et al. 2006; 
Obergaulinger et al. 2006; Bisnovatyi-Kogan, Moiseenko & 
Ardelyan 2008; Benson & Babul 2009). 

Processes that lead to in-situ amplification of large-scale 
magnetic field, due to the dynamics of the flow, are usually 
referred as dynamos (Parker 1955, 1987; Moffatt 1978; Bur- 
bidge, Layzer & Phillips 1981; Zeldovich & Ruzmaikin 1987; 
Roberts & Soward 1992; Kulsrud et al. 1997; Sato 1999; 
Brandenburg & Subramanian 2005). There are many kind 
of dynamo processes but, essentially, all involve twisting of 
the magnetic fieldlines by the flow and reconnection events 
that allow to rearrange irreversibly the field topology. Dy- 
namos usually require a fully three-dimensional flow struc- 
ture. Cowling's theorem, for example, states that dynamos 
cannot work for two-dimensional flow patterns. 

Of particular interest in astrophysics are dynamo pro- 
cesses due to the presence of small-scale fluctuations in the 
flow and turbulence. The idea that small-scale velocity and 
magnetic field fluctuations might be correlated, leading to a 
large-scale effective electromotive force capable of amplify- 
ing and generating a large-scale magnetic field, is at the base 
of the so-called mean-field dynamo theory (Moffatt 1978; 
Krause & Raedler 1980; Brandenburg & Dobler 2002; Bran- 
denburg & Subramanian 2005). Mean-field dynamos have 
been applied to a large variety of astrophysical systems, from 
the Sun (Parker 2009) to stellar magnetism (Brandenburg 
& Dobler 2002), from accretion disks (Khanna & Camen- 
zind 1996b; Pariev, Colgate & Finn 2007) to proto-neutron 
stars (Bonanno, Rezzolla & Urpin 2003), from the origin of 
the galactic field (Shukurov 2002) to that of the cosmologi- 
cal primordial field (Kulsrud & Zweibel 2008), just to cite a 
few. 

The formulation of a mean-field closure of Maxwell's 
equation in the context of General Relativity (GR) has been 
done only for two astrophysical cases, to our knowledge: 
Marklund & Clarkson (2005) have investigated the origin of 
the cosmic magnetic field during inflation, while Khanna & 
Camenzind (1996a) and Khanna (1999) have focused on the 
specific case of disks in Kerr metric. Both have performed 
some study in the limiting kinematic case, where the back- 
reaction of the magnetic field on the plasma is neglected. 
Khanna (1999) also have not considered the displacement 
current, which however might become non-negligible near 
the event horizon of a BH, or for rapidly evolving systems. 



They show that various new terms can arise in GR which 
are not present in a flat space-time. There is also a debate 
if Cowling's anti-dynamo theorem still holds in GR: frame 
dragging effects can generate currents even for axisymmetric 
configurations (the gravito-magnetic term) in the absence of 
turbulence (Khanna & Camenzind 1996a). 

The possibility of a stable and reliable mean-field clo- 
sure for turbulent dynamos is however very important in 
the context of numerical simulations. Quite often large- 
scale simulations are required to investigate the dynamics of 
GRB engines, accretion disks in AGN, and magnetosphere- 
jet coupling. The development of small-scale fluctuations is 
a property of turbulence that make its numerical investiga- 
tion within global models prohibitive if not impossible at 
the moment. The idea behind a mean-field approach is that 
the effects of physical processes at scales that cannot be re- 
solved can instead be modeled by an appropriate closure 
of the equations, so that the problem can become treatable 
by numerical investigation. However we want to stress here 
that, in the mean-field approach, determining what could 
be realistic values for the parameters that are used in the 
closure, is usually non trivial, and often requires the use of 
mesoscale informations, and extrapolation of flow proper- 
ties to small unresolved scales. It rests to be proved that 
a mean field approach can achieve full resolution of micro- 
scopic physics on the macro-scale. 

Here we present the first fully covariant mean-field dy- 
namo closure of Maxwell's equations, by extending the co- 
variant Ohm's law widely used in resistive GRMHD (Lich- 
nerowicz 1967; Anile 1989). The new contribution is due to 
an a-dynamo term proportional to the large-scale magnetic 
field as measured in the local comoving frame, in analogy to 
the classical case. Moreover, having in mind numerical ap- 
plications, the equations are then cast in the 3 + 1 formalism, 
leading to a modified evolution equation for the spatial elec- 
tric field as measured by Eulerian observers. When dynamo 
effects are negligible, the equation reduces to that already 
derived for resistive MHD in special relativity (Komissarov 
2007; Palenzuela et al. 2009; Dumbser & Zanotti 2009), 
thus extending it to the more general case of 3 + 1 resistive 
GRMHD. It is not our intention to discuss here either the 
validity of mean-field dynamo or its theoretical implication 
within GRMHD. As already done for the purely resistive 
case, in the above cited works, here we mainly focus on the 
numerical implementation of the proposed closure. Our nu- 
merical references are the ECHO and X-ECHO codes (Del 
Zanna et al. 2007; Bucciantini & Del Zanna 2011), and the 
actual implementation and validation tests of this new dy- 
namo closure will refer to these schemes. 

The plan of the paper is the following. In Sect. 2 we 
discuss the mean-field closure for the induction equation, 
and present our covariant GR extension. Sect 3 is devoted 
to its representation in the 3 + 1 formalism, necessary for 
numerical modeling, whereas in Sect. 4 we discuss its actual 
implementation within the X-ECHO code for GRMHD. In 
Sect. 5 we present a set of simple standard tests done both in 
the so called kinematic and in the fully dynamic regime, and 
compare them with previously published results and with 
analytic and semi-analytic solutions. Finally we present our 
conclusion in Sect. 6. 

In the following we assume a signature —,+,+,+ for 
the space-time metric and we use Greek letters fx, u, A, ... 
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(running from to 3) for 4D space-time tensor components, 
while Latin letters i,j, k, ... (running from 1 to 3) will be em- 
ployed for 3D spatial tensor components, and spatial vectors 
will be often written using bold face characters. Moreover, 
we set c = G — Mq = 1 and all factors will be absorbed 
in the definition of the electromagnetic fields. 



2 THE COVARIANT OHM'S LAW AND ITS 
MEAN-FIELD DYNAMO CLOSURE 

Let us now briefly discuss the main idea behind the mean- 
field dynamo theory (Krause & Raedler 1980). Ohm's law 
for resistive (classical) MHD reads 

E + vxB = r,J; J = VxB, (1) 

where, E, B, v, J are the electric field, the magnetic field, 
the velocity and the current respectively, and rj is the resis- 
tivity or coefficient of (isotropic) magnetic diffusivity. If now 
one separates the quantities into their large-scale mean val- 
ues E, v, B, and small-scale fluctuating parts SE,Sv,5B, 
a new electromotive force due to turbulent motion appears 
in Ohm's law 

E + v x B = -Sv x SB + r/ J; J = V x B. (2) 

In general the small-scale fluctuating quantities are corre- 
lated and thus their product has a non- vanishing mean. The 
key assumption is that this mean can be written as a func- 
tion of the mean quantities and their derivatives, and in the 
simplest case that it is linear in the value of both the mean 
magnetic field and its curl. Namely, it is often assumed that 

Sv x SB = aB — /3V x B (3) 

where the two scalar coefficients are both proportional to 
the local turbulent correlation time r c . In particular 

a = -\t c 5v X (V X Sv); /3=±t c 5^, (4) 

where the a-term is proportional to the kinetic helicity and 
the /3-term is related to a random walk for fluid elements. 
For a deeper insight on the properties of turbulence and a 
more exhaustive derivation of the a and /3 terms the reader 
is referred, for instance, to Kulsrud (2005). Inclusion of mag- 
netic helicity in the definition of a and a different mean-field 
dynamo closure allowing for a dynamical definition of the 
electromotive force due to small-scale turbulent fluctuations 
may be found in Blackman & Field (2002). 

Dropping the bars and referring from now on to just 
large-scale averaged quantities, the classical form for the dy- 
namo closure is then 

E + v x B — —a B + (/3 + rj) J . (5) 

Notice that in general both a and /3 will be tensors, how- 
ever we will focus here on the isotropic case where they 
can be dealt with as scalars. Their values might depend on 
fluid quantities like density, temperature, or magnetic field 
strength. Moreover, even the specific physical problem at 
hand could have an influence, as the resulting asymptotic 
turbulent state may be strongly affected by the assumed 
initial conditions. When these coefficients can be treated 
as constants, the induction equation for classical MHD be- 
comes 

dtB = Vx(i)xB) + nVxB + (/J + tj) V 2 B, (6) 



and it is now apparent that the presence of the a-term may 
introduce exponentially growing modes that are known as 
mean-field dynamo waves, whereas the /3-term acts as a sort 
of turbulent diffusivity or turbulent resistivity, which is often 
dominant over the kinetic one (in the fast dynamo case). 
The /3-term can be interpreted as due to turbulent mixing: 
the convection cells mix up magnetic field lines of different 
polarities on small scales, and thus reduce the mean field, 
which is equal to the field averaged over larger scales In the 
kinematic regime, a, /3 (and rj) are input parameters, as well 
as v, and only the above equation needs to be solved. 

Let us finally summarize mean-field dynamo treatment 
within classical MHD by rewriting the classical Ohm's law 
above as 

E'=ZB + r)J. (7) 

that is by replacing E + v x B — > E' , the electric field in the 
frame comoving with the fluid, — a — > £ (to avoid conflict 
with the lapse function in the 3 + 1 standard notation, to be 
introduced in next section), and /3 + n — > 77, combining mag- 
netic and turbulent diffusivity in a single coefficient. In the 
remainder of this paragraph we will propose a fully covari- 
ant generalization of Eq. 7, which is novel in the literature, 
to our knowledge. 

The covariant Maxwell's equations are written in terms 
of the Faraday (antisymmetric) tensor F^ v and its dual 
F*^ = l e ^ XK F XK , where e^ XK is the spacetime Levi-Civita 
pseudo-tensor (here we use the convention (— g) 1 / 2 e 0123 = 
-(-ff)~ 1/2 eoi23 = 1), as 

V M F* Ml - = 0, V M F^ = -/", (8) 

where i 7 * is the 4-current. The above quantities may be de- 
composed in the reference frame comoving with the fluid 
4- velocity u M as 

F"" = u l "e v -e^u" + €^ XK bxu K , (9) 

F" 4 " = un v - - e^ XK e x u K , (10) 

and 

I» = q u»+f, (11) 

where e" = F*"u v , fe" = F*^ v u v , q = -J"u M , and f are, 
respectively, the electric field, magnetic field, charge density, 
and (conduction) current measured in such frame (e M tt M = 
= = 0). 

Ohm's law for (isotropic) resistive GRMHD is usually 
written as a linear relation between the comoving electric 
field and current (Lichnerowicz 1967; Anile 1989) 

ef l = r,f, (12) 

and the ideal GRMHD relation of a vanishing comoving elec- 
tric field e M = is recovered by letting r\ = (an ideal 
plasma with infinite conductivity), which was the closure 
employed in Del Zanna et al. (2007) and Bucciantini & Del 
Zanna (2011). The straightforward extension to include a 
mean-field a-dynamo effect appears to be 

e" =( > b lx + r 1 f, (13) 

where a new term proportional to the comoving mean mag- 
netic field 6 M now appears and again only the isotropic case 
has been considered. Both the coefficients £ and r\ serve as a 
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sort of sub-grid modeling of the turbulent motions, and tur- 
bulent diffusivity is supposed to be in general higher than 
its kinetic value, as in the classical case. Eq. 13 is thus our 
proposed fully covariant generalization of Eq. 7, to which it 
correctly reduces in the comoving frame where v — (see 
next section). 

We want to stress here that, as in the standard mean- 
field dynamo theory, there is a large degree of freedom in the 
choice of the closure relations. It can even be debated if a 
closure in terms of mean quantities and their derivative can 
be found at all. The quest for an appropriate closure is only 
apparently more complex in general relativity, due to the 
requirement of general covariance, but practically this seems 
to support the validity of our simple expression in Eq. 13. 
However, we want to stress once again that this can be so 
straightforward only when the isotropic case is assumed, as 
done in this work for simplicity. Anisotropic resistive MHD 
in a Minkowskian spacetime was considered by Zanotti & 
Dumbser (2011). On the other hand, our scalar parameters 
£ and rj may be both function of all the other (macroscopic) 
quantities, and thus evolve dynamically in time with them. 
For simplicity, however, even if the scheme is built to take 
into account this most general case, numerical tests will be 
presented only for constant £ and rj. 

3 THE DYNAMO CLOSURE IN 3 + 1 GRMHD 

In the 3 + 1 formalism the line element is usually given as 

ds 2 = -a 2 dt 2 + ^ 3 {dx l + P l dt){dx 3 + /3 j dt), (14) 

where a the lapse function, f3 % the spatial shift vector, both 
arbitrary due to gauge invariance in the choice of coordi- 
nates, and 7ij is the spatial 3-metric, with determinant 7 
(— g — a 2 -y). In this metric the unit vector of the Eulerian 
observer's 4-velocity is 

n M = (-a, 0), n" = (1/a, -/37a), (15) 

and projection onto the spatial hyper-surfaces normal to n M 
is achieved via the 3-metric 

Ivv = g» v +n^n„. (16) 

Spatial projection for a generic vector V M (or tensor) is then 
_L V = yZV = V M + {y v n v )n' i , and for a spatial vector 
_L V M = . Such a spatial vector must have a vanishing 
contravariant temporal component, since V^n^ — 0. 

If we now decompose the electromagnetic quantities 
within the 3+1 Eulerian framework, in analogy with the 
previous relation we write 

F"" = n"B" - fi"n" + e^ XK B x n K , (17) 

F*»" = n»B v - B»n v - e^ XK E x n K , (18) 

and 

I M = gn M + J M , (19) 

where = F M "n„, B M = F*^n v , q, J M are, respectively, 
the electric field, magnetic field, charge density, and (con- 
duction) current measured in such frame (i? M n M = B^n^ = 
J^n^ — 0). The Maxwell's equations take the usual form, 



plus some extra terms due to 3 + 1 GR metric 

■y- 1/2 dt(j 1/2 E)-V x (qB-/3x E) = -(aJ -qp), (20) 

7 _1/2 9i(7 1/2 B) + V x (aE + P x B) = Q, (21) 

V-E = q, (22) 

V • B = 0, (23) 

and we do not repeat here the momentum-energy conserva- 
tion equation, the continuity equation for the mass density 
(an equivalent one holds for the electric charge density q), 
and Einstein's field equations in the 3 + 1 formalism. Com- 
pared to ideal GRMHD, where the electric field is a derived 
quantity, we must now also consider the evolution and con- 
straint equations for E, where the sources q and J appear. 

Let us see now how to treat the various forms of the 
generalized Ohm's law. It is first convenient to decompose 
the quantities related to the frame comoving with the fluid 
within the 3+1 Eulerian split of time and space, namely 

= Fn" + Fv" , (24) 

e M ^T{E-v)n" +Y{E'" +e" vX v u B x ), (25) 

b" = T(B ■ v)n" + T(B M - ^ vX v v E x ), (26) 

f = (q- g rK + r - qoTv*, (27) 

where T = (1 — v 2 )^ 1 ^ 2 is the Lorentz factor of the fluid 
flow, go = T(q — J ■ v), and e'" /X = e^ vXn n K is the spatial 
Levi-Civita pseudo-tensor, for which 7 1 / 2 e 123 = 1. We now 
have three possibilities: 

(i) Ideal GRMHD (e M = 0): 

The spatial projection readily provides the usual ideal 
MHD assumption 

E + v x B = 0, (28) 

exactly as in the classical limit. 

(ii) Resistive GRMHD (e" = J?j M ): 

The time projection leads to F(E ■ v) — r/(q — qoF), which 
may be used to express qo in the spatial component. The 
result is 

T[E + v x B - (E- v)v] = rj(J - qv), (29) 

as already found in previous treatments within special rel- 
ativists resistive MHD (Komissarov 2007; Palenzuela et al. 
2009; Dumbser & Zanotti 2009), thus here we extend its 
validity to 3 + 1 GRMHD. When ij = 0we reduce to the 
previous ideal MHD case, while for |u| <§C 1, \E\ ~ M|-B| wc 
recover the classical limit of Eq. 1. 

(iii) Resistive GRMHD + dynamo (e M = £6 M + ??j M ): 
The time projection is now 

r(E.v) = r](q-qor) + ^(B-v), (30) 

which again may be used to express go in the spatial com- 
ponent. Now the result is 

T[E+vxB-(E-v)v] = r)(J-qv)+££[B-vxE-(B-v)v\, 

(31) 

which is a novel closure to our knowledge, reducing to the 
case of resistive GRMHD when £ = and to the classical 
case of Eq. 7 for small velocities. 

It is interesting to note that in all cases the presence of 
a curved or even evolving GR metric is not apparent in 
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Maxwell's equations, since a or /3 1 terms do not appear ex- 
plicitly, whereas 7^ is just needed to work out scalar and 
cross products between (spatial) vectors. 

While in resistive schemes for classical MHD [see e.g. 
Landi et al. (2008)] q and E do not play a role and the 
system is closed simply by taking J — V x B in Ampere's 
law, the presence of Maxwell's displacement current in the 
relativistic case forces one to use Eq. 20 to evolve the electric 
field in time. Only in ideal GRMHD is it still possible to 
neglect the Maxwell equations for the electric field, since E 
is provided from the ideal Ohm's law Eq. 28. In the other 
cases the two equations for E must be evolved or preserved 
in time, and we need to face the problem of dealing with the 
unknown sources q and J. In the proposed scheme we use 
Ohm's law to provide an expression for the spatial current 
density J, while the constraint on q is enforced by using the 
Gauss theorem in Eq. 22. Notice that here we do not evolve 
the charge density via the corresponding continuity equation 
(charge conservation law), since we choose to impose directly 
q = V • E. 

In the most general case, including mean-field dynamo 
effects, the result is 

^ 1/2 d t (-y 1/2 E) - V x (aB -/3xE) + (av- (3)q = 
- aV[E + vx B (E- v)v]/r] 

+ £aT[B -vx E - (B ■ v)v]/r). (32) 

When 77 = we can neglect the time evolution term and 
the rest of the left hand side, therefore time integration of 
the electric field is not required; when also £ = we recover 
the ideal Ohm's law condition E = —v x B, as expected. 
When r\ > problems arise because Eq. 32 is usually stiff, 
especially for low values of the resistivity, since the source 
terms on the right are larger than the time marching term 
by a large factor I/77, and some sort of implicit numerical 
scheme must be employed for time integration. Note that 
in this respect the presence of an dynamo a— term adds no 
further complexity to the numerical algorithm, and its im- 
plementation within an existing resistive code is expected to 
be straightforward. 

4 IMPLEMENTATION IN THE ECHO CODE 

Let us discuss here how the closure relation Eq. 32 can be 
solved in a numerical scheme for 3 + 1 GRMHD taking as a 
reference the ECHO (Del Zanna et al. 2007) and X-ECHO 
(Bucciantini & Del Zanna 2011) codes developed in the ideal 
regime. Here we propose a very simple method to integrate 
the equation implicitly. Considering a first order discretiza- 
tion of the time derivative, one can write an expression for 
the electric field at the end of a timestep At as a function 
of other quantities, both at the end (superscript ^) and at 
the beginning (superscript ^) of the implicit procedure. 

Introducing the spatial components instead of the vec- 
torial notation we have 

E iW =7 (i)- 1/2 7 (0) 1 / 2 Q*(0)_ 

[r (D^(D + e «* e (i) B (U _ (E k ^v k 1) )v iW /r^]/fj+ 

ZF^BW - e ijk v^E^ - (B k ^v^)v^/r^]/fj (33) 

where rigorously fj = r)^ and £ = since the two coeffi- 
cients might be in principle function of other variables like 



temperature, density or magnetic field. The other assump- 
tions are 

v — IV, Vi = Tvi, r 2 = 1 + ViV 1 , 

Q* = E l + At[-(av l - p)q + ^ k d,(aB k - e klm /3 l E m )\, 

(34) 

q = ^ 2 d k (^ 2 E k ), 
1/fj = At a/rj. 

If we now solve for E l , after some lengthy algebra we find 
E i [T + ^ + £ 2 (r 2 -l)/(r + 77)] = 

- t ij %B k + fj[Q' + {Q k v k )v l ]/{1 + fjT) 
+ £[RB l -fj(B k v k )v i /(l+fjT)] 

- £[(r 2 - 1)B> - {v k B k )v l + fje^vjQkyiT + fj) 

+ eir^ k v J B k ]/(r + f l ) 

+ e[rir(Q k v k ) + aB k v k )]v/[(i + f,r)(r + fj)], (35) 

where for simplicity we have let 7 (1) 1/2 -y( ) 1/,2 Q l («) _> g» 
and dropped all superscripts. When £ = 0, that is in the 
purely resistive GRMHD case, many terms cancel out and 
we are left with the simple relation 

E\V + fj) = -e ijk VjB k + f,[Q l + (Q k v k )v l ]/(1 + fjT), (36) 

which automatically includes the limit for an ideal plasma 
E i = -e iik VjB k when r) = 0. 

4.1 Primitive variables 

Eq. 35 provides the electric field at the end of a timestep as 
a function of other primitive variables like the velocity and 
the magnetic field. However numerical schemes for GRMHD 
usually evolve conserved variables like momentum and en- 
ergy, and primitive variables are not immediately available 
at the end of a timestep, but must be derived by inverting 
a set of non-linear equations. This implies that Eq. 35 must 
be solved simultaneously with the inversion from conserved 
variables to primitive. In analogy to Palenzuela et al. (2009), 
the derivation of the primitive variables and the electric field 
is done along the following lines: 

• Given the conserved variables at the end of a timestep, 
a guess for the pressure p* is chosen, using the value at the 
previous timestep. 

• With this value p* of the pressure kept fixed, a guess 
for the v l * is chosen, again using the values at the previous 
timestep. 

• E % components are derived according to Eq. 35. This 
step is performed only when the solution of an implicit equa- 
tion is required. In general the Q % contain all of the explicit 
terms (see the following subsection on time-stepping) . Eq. 34 
provides the value of Q l in the simple case for a first order 
implicit solver. 

• The momentum equations are inverted keeping fixed 
the value p* , by means of a Newton- Raphson scheme, where 
the Jacobian is computed numerically, to provide a new 
guess v % * for the v* components. A new electric field is de- 
rived and this loop is iterated until convergence. 

• The energy equation is finally inverted by means of a 
Newton-Raphson scheme using the values of the velocity v % * 
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and electric field obtained in the inner cycle, to provide a 
new guess p* for the pressure. Again the derivative required 
for the Newton-Raphson scheme is computed numerically, 
allowing for a general EoS. 

• The overall cycle over the pressure p is repeated until 
convergence. 

This approach has the advantage to allow the use of 
a general equation of state, hence it is not limited to the 
ideal gas or polytropic EoS, and in principle can be general- 
ized to other closure relations for Ohm's law. We opted for 
numerical Jacobians, as opposed to approximations to the 
true analytical ones (Palenzuela et al. 2009), which might 
be extremely complex to derive and expensive to compute. 

The above implementation is straightforward in the 
Cowling approximation, where the metric terms are fixed 
in time, but it can be easily extended also to a dynami- 
cal space-time by using for example the XCFC approach as 
used in X-ECHO (Cordero- Carrion et al. 2009; Bucciantini 
& Del Zanna 2011), given that the equation for the evolu- 
tion of the metric terms are decoupled from the inversion 
from conserved to primitive variables. The only problem is 
actually the treatment of the lapse function, appearing in 
the definition of fj. While in XCFC the conformal factor and 
the determinant of the three-metric 7 can be easily com- 
puted from the conserved quantities, the lapse a requires 
the previous knowledge of the primitive variables, including 
the electric field. This implies that one should in principle 
solve simultaneously for a and the primitive variables. Given 
that in general the lapse is a slowly and smoothly varying 
function of time, we prefer to use for simplicity a' ' in 77, 
and we expect this choice to introduce only minor errors. 



4.2 Time stepping and constrained transport 

The scope of this work is to present and verify the imple- 
mentation of a dynamo closure in a numerical scheme for 
3 + 1 GRMHD. Our code allows the use of two distinct time- 
stepping approaches: a simple lst-2nd splitting of the tem- 
poral evolution, between the implicit and the explicit part, 
and a more rigorous 2nd order IMEX Scheme (Pareschi & 
Russo 2005; Palenzuela et al. 2009). Let us describe here 
both of them with reference to a typical system of stiff and 
non-stiff equations, for the conserved variables U, W: 

d t U = F(U,W), 

d t W = G(U,W) + ^R(U,W). (37) 

In the simple lst-2nd scheme, the explicit part is solved 
using a modified Eulerian scheme which is second order in 
time, while a simple first-order scheme is used to solve the 
implicit one. 

jjW = jjn + AtF(U",W n )/2, 

w m = W " + AtG{U",W n )/2 + —R(U W ,W W ), 
U n+1 =U n + AtF(U m ,W {1) ) (38) 
w n+i = w n + AtG(t7 (i) ) W W) + ^R(U n+1 ,W n+1 ), 

where the implicit step acts only on the electric field E, 



it is peformed during the inversion from conserved to primi- 
tive variables (see above), accoding to the solution provided 
in the previous section. 

This symple scheme has the advantage of being easily 
implementable, requiring only a modification in the defini- 
tion of the electric field, within the algorithm that derives 
the primitive from the conserved variables. Moreover the al- 
gorithm is well behaved in the case 77 = [it reduces to 
solving R(U, W) = 0], and thus it can handle also the Ideal 
MHD regime. 

The second order IMEX that we have implemented is 
the following: 

w m = w n + ^tt R{uW wW) 
v 

U {2) =U n + AtF(U w ,W w ), 
W (2) = w ™ + AtG(U (1) ,W w )+ 

+ y [(1 - 2f,)R(U w ,W w ) + nR(U (2 \ W {2) )\, 

jjn+i = u n + ^* [F ( t7 (i) > W (D) + F (U (2 \W {2) )}, 

w n+l =W n + ^ [G([7 (l) j W W) + G (U {2) ,W {2) )} + 

+ ^ [R{U W )W W ) + R(U {2 \W (2) )}, (39) 
where /u = 1 — l/y/2. 

Again the implicit step acts only on the electric field 
E, it is peformed during the inversion from conserved to 
primitive variables. Note however that now the Q % are no 
more given simply by Eq. 34. 

However, in comparison with the previous lst-2nd 
scheme, this IMEX scheme requires 3 different steps instead 
of 2, the steps do not have the same functional form, and, 
because of the last step, it is not well behaved in the case 
77 = 0, and so it cannot be used for Ideal MHD problems. 

All the test that we present in the following have been 
repeated both with out lst-2nd scheme and with the IMEX 
scheme. As we will discuss, depending on the problem or the 
desired accuracy, our proposed lst-2nd scheme might offer 
and easy alternative to more sophisticated IMEX implemen- 
tations. 

The value of the timestep At is chosen in accordance to 
the CFL condition, with typical Courant numbers ranging 
from 0.1 to 0.5 (smaller values have been adopted in high re- 
sistivity runs to avoid large truncation errors in the solution 
of the implicit part, due to the smaller diffusive timescale). 

As far as the electromagnetic constraints of Eqs. 22-23 
are concerned, all previously developed schemes for (special) 
relativistic resistive MHD (Komissarov 2007; Palenzuela 
et al. 2009; Dumbser & Zanotti 2009) adopt a divergence- 
cleaning approach (Munz et al. 2000; Dedner et al. 2002), 
where an augmented system of equations is introduced 
(including the charge conservation law), and where the 
solenoidal constrain on the magnetic field and Gauss's theo- 
rem are not enforced but preserved by damping and propa- 
gating away any violation arising from numerical truncation 
errors. Here instead we choose to use a fully constrained 
scheme, where the charge q appearing in the equation for 
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E is taken directly from Eq. 22 while the divergence-free 
condition on B is treated with staggered grids via the Up- 
wind Constrained Transport method (Londrillo & Del Zanna 
2000; Londrillo & del Zanna 2004; Del Zanna, Bucciantini 
& Londrillo 2003a). A benefit of this fully constrained ap- 
proach is that, neglecting dynamo effects (£ = 0) in the 
purely resistive case, it is possible to recover the ideal MHD 
expression simply by setting n — (see Eq. 35), while all 
previous scheme recover the ideal regime only as a limit for 
low resistivity r\ — > 0. However this also depends on the 
time-stepping algorith, as we have discussed previously. 

To conclude this section, note that in the resistive case 
the maximum wave speed is not limited by the fast mag- 
netosonic mode and may approach the speed of light inde- 
pendently of the value of the magnetic field. The Riemann 
problem at cell interfaces is thus solved for simplicity us- 
ing the maximally diffusive (global) Lax-Friedrichs scheme, 
where the fastest characteristic speed is set to be equal to the 
speed of light. The use of such a scheme, as opposed to more 
accurate Riemann solvers like the HLL or HLLD, might lead 
to less satisfactory results in the ideal MHD limit, where 
sharp discontinuities are allow to arise. In the resistive case, 
where smooth profiles are expected, the use of more diffu- 
sive algorithms is less problematic, especially in conjunction 
with high-order reconstruction algorithms, a distinguishing 
feature of our ECHO scheme. Again, we leave the implemen- 
tation of more accurate Riemann solvers as a future upgrade. 



5 TEST PROBLEMS 

In the following we present a set of standard tests. The first 
three are done in the purely resistive regime, for compari- 
son with results previously presented in the literature. Then 
there are three dynamo problems in the so-called kinematic 
regime, and finally a fully dynamical problem All but the 
last of these tests are done in a flat, stationary metric, since 
they are aimed at evaluating the implementation of the re- 
sistive/dynamo closure, and for a more straightforward com- 
parison with previous results, and with analytical solutions. 



5.1 Resistive tests 

We present here a set of three resistive tests (£ = 0) in a 
stationary Cartesian grid (Minkowskian spacetime), to be 
compared with both analytical and previously published re- 
sults. In the following we use a third-order CENO spatial 
reconstruction with MC limiter, and both our fst-order im- 
plicit / 2nd-order explicit temporal evolution scheme, and 
the 2nd order IMEX scheme. A maximally diffusive global 
Lax-Friedrichs Riemann solver (LF) is employed for upwind- 
ing. 



5.1.1 Self-similar current sheet 

This problem was first proposed by Komissarov (2007), and 
it has been presented also by Palenzuela et al. (2009) and 
Dumbser & Zanotti (2009). It is a truly resistive problem, 
which allows for the analytical self-similar (in t/x 2 ) solution 



Reslsitive CS 




X 

Figure 1. Evolution of the magnetic field in a self-similar current 
sheet, for r\ = 0.01. The red dot-dashed line is the initial condition 
at t = 1. The blue solid line is the numerical solution at t = 
10, indistinguishable from the green dashed line representing the 
exact solution Eq. 40. 

in the limit of infinite pressure: 

^ (M )=iw(-y (40) 

where erf is the error function. 

Despite the evolution being almost purely resistive (in 
principle, in the limit of infinite pressure only the mag- 
netic field evolves and only the induction equation needs 
to be solved), the problem is followed in the fully dynamical 
regime. The initial conditions at t = 1 are: p = 1, p = 50, 
E i = 0, v l = 0, B x = B z = 0, B y = B y (x, 1), with B = 1, 
and we have adopted an adiabatic coefficient 7 = 4/3. Our 
computational domain extends in the range x — [—1.5, 1.5], 
and is covered by a uniform grid with 200 cells. The prob- 
lem is evolved to the final time t = 10. In Fig. 1 we compare 
our numerical results with Eq. 40, in the case n — 0.01. 
Errors and convergence estimates are presented in Sect. 5.3 
together with a comparison between the lst/2nd scheme and 
the IMEX scheme. 

5.1.2 Resistive shock tube 

This problem was proposed by Dumbser & Zanotti (2009), 
and differs from the one presented in Palenzuela et al. 
(2009), which is done in the transverse MHD regime. Shock 
tubes are excellent tests for monitoring the shock-capturing 
properties of a numerical scheme. The initial conditions are: 

(p,p,v x ,v y ,v z ,B x 1 B y ,B z ) = 

(1.08, 0.95, 0.4, 0.3, 0.2, 2.0, 0.3, 0.3) for x < 0. (41) 

and 

(p,p,v',v",v* i B",B",B*) = 

(1.0, 1.0, -0.45, -0.2, 0.2, 2.0, -0.7, 0.5) for x > 0. (42) 

and the problem is followed to a final time t = 0.55. The 
initial electric field is set equal to the Ideal MHD value 
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Figure 2. Resistive shock tube problem. The upper panel shows 
the density and the lower panel shows the y-component of the 
magnetic field at the final time t = 0.55. The blue solid line is the 
case rj = 0, the dotted cyan line is the case r\ = 0.01, the green 
dashed line is for rj = 0.1, the yellow dot-dashed line for rj = 1, 
while the red dash-triple dotted line refers to r\ = 1000. 



E' — —e^ BjVk. The computational domain extends in the 
range x = [—0.5, 0.5] and is covered by a uniform grid with 
400 cells. The test was repeated with the following values 
for the resistivity: r\ — 0,0.01,0.1,1,1000, where the last 
value was chosen to be big enough such that the evolution 
practically corresponds to the zero conductivity case. The 
adiabatic coefficient is here 7 = 5/3. Results are shown in 
Fig. 2. 



5.1.3 Resistive rotor 

This is a fully multidimensional problem. The relativistic 
Ideal-MHD version was first proposed by Del Zanna, Buc- 
ciantini & Londrillo (2003b), while the resistive version has 
been presented by Dumbser & Zanotti (2009). 

A circular region with radius r = 0.1, uniform den- 
sity p = 10, and rotating with a uniform angular velocity 
$7 = 8.5 is located within a medium at rest, with a lower 
density p = 1. The pressure p = 1 and the magnetic field 
(B x , B y , B z ) — (1,0,0) are uniform in the whole domain. 
The adiabatic coefficient is 7 = 4/3 and the system is fol- 
lowed to a final time t — 0.3. The initial electric field is set 
equal to the ideal MHD value E 1 = ~e ljk BjV k . The compu- 



tational domain is x — [—0.5,0.5], y = [—0.5,0.5], with the 
rotating region located at the center, and is covered with 
a uniform grid of 400 x 400 cells. The problem is solved 
in the ideal MHD regime rj — 0, in the quasi-ideal regime 
rj = 0.001, and in the resistive regime rj — 0.1, and the 
results are shown in Fig. 3. Comparison between the case 
rj = and rj — 0.001 suggests that, at this resolution, the 
intrinsic resistivity of the scheme is smaller than rj = 0.001. 
Moreover the case rj = 0, when repeated with a HLL solver 
(not shown here) does not show any significant improvement 
with respect to the same case done with LF. 



5.2 Dynamo tests 

Here we present the tests done using the closure Eq. 13, lead- 
ing to the full version in Eq. 32 of the evolution equation for 
E. The first three are performed in the so-called kinematic 
regime, where only the induction equation is solved and only 
the electric and magnetic fields are allowed to change in 
time. Given the physical nature of mean-field dynamos, re- 
lated to small-scale turbulence, the kinematic approximation 
is actually commonly adopted as a suitable approximation. 
In fact, the magnetic field generated by turbulent dynamo 
action is often well below equipartition, with respect to the 
internal energy density, and as such of negligible dynamical 
consequences. Moreover, kinematic dynamos often allow for 
simple analytical solutions, whereas the non-linear feedback 
on the flow can only be dealt with using numerical methods. 
As for the resistive cases, in all of the following we employ 
a third-order CENO spatial reconstruction with MC lim- 
iter, the maximally diffusive Lax-Friedrichs solver, and both 
the usual lst-order implicit / 2nd-order explicit and IMEX 
time-stepping schemes. 



5.2.1 ID steady dynamo 

This is a very simple test describing the growth of magnetic 
field in a stationary medium. It is the dynamo equivalent of 
the resistive current sheet presented in section 5.1.1, because 
it admits an analytical solution (see Appendix A). 

We perform a few runs with different resistivity rj = 
0.05,0.1,0.25, different values of the wave number k, while 
£ = 0.5 is kept fixed. The computational domain is x = 
[— tv, 7r] is covered with 200 uniformly spaced zones, and the 
problem is followed to a final time t = 100. Our initial con- 
ditions are chosen to correspond to the growing eigenmodes 
of the problem (see Appendix A): 



B v = O.lsin(fcr), B x = 0, 
B z = -0.1 cos (fcc) 



(43) 



Fig. 4 shows the evolution of the B y component of the 
magnetic field, compared with the corresponding analytical 
expectations, for various values of rj and k. It is interesting 
to note that in the case rj = 0.1 and k — 5, the magnetic 
field is expected to remain unchanged, due to the opposite 
effects of resistivity and dynamo. Given the nature of the 
dynamo closure in Eq. 13, due to round-off and interpolation 
errors, the evolution of the fastest growing mode however 
dominates asymptotically. This is what is observed around 
t = 60: the field starts to grow exponentially, with a growth 
rate corresponding to the fastest growing mode. The exact 
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Figure 3. Relativistic rotor. The upper panels show the pressure, while the lower panels show the z-component of the electric field at 
the final time t = 0.3. Left column: the Ideal case r\ = 0. Central column: quasi-ideal case r\ = 0.001. Right column: resistive case w = 0.1. 



time, when this transition happens, depends on roundoff and 
interpolation errors, and differs for the IMEX scheme. 

In Fig. 5 this property of the dynamo solution is shown 
for 77 = 0.05: the initial solution with a wave number k = 9 
evolves with a slow growth rate until the fastest growing 
mode, corresponding to k = 5, emerges at t — 12, to dom- 
inate the subsequent evolution. Note that, given the expo- 
nential nature of the solution, the transition is very sharp. 
Errors and convergence estimates are presented in Sect. 5.3 
together with a comparison between the 1st /2nd scheme and 
the IMEX scheme. 



5.2.2 Thin shear layer 

In this section we present a two-dimensional shear dynamo 
problem in the so-called thin-layer approximation, where one 
of the dimension is assumed to be much smaller than the 
other, such that higher order variations of all quantities in 
that direction can be neglected. This is the relativistic equiv- 
alent of the ID problem investigated by Arlt & Riidiger 
(1999), where spatial variations of the dynamo coefficient 
and quenching were also included. It is immediately evident 
the difference between the non-relativistic and the relativis- 
tic case. In the former, one just need to add a shear term 
to the induction equation, proportional to the derivative of 
the shear velocity along the neglected dimension. However, 



this is not possible within our relativistic formalism, where 
the modification due to the dynamo process does not appear 
in the induction equation, which for us retains the general 
form given by Maxwell's equations, but instead in the defi- 
nition of the electric field via the closure of Eq. 13. For this 
reason, we adopt the thin-layer approach that allows one to 
recover the ID results in the limit of a negligible thickness. 
This test is representative of typical conditions in accretion 
disks, in the limit where one models the differential rotation 
with a shear term in the radial direction and consider modes 
extending along the vertical direction. 

The computational domain is x = [ — 7r, 7r], with periodic 
boundary conditions, and z — [—5,8], with boundary con- 
ditions where all quantities are linearly extrapolated. The 
resistivity is set 77 = 0.1, while the dynamo term is £ = 0.2. 
A shear velocity is imposed to be v y — 0.9z. The problem is 
followed to a final time t = 24, corresponding to a dynamo 
wave period. In Appendix A the relativistic solution is de- 
rived, and we provide also the numerical results for the set 
of parameters of the present test. Our initial conditions are 
chosen to correspond to the expected growing eigenmode of 
the problem 

B v = 0.212 cos(a; - 0.662), B x = 0, B z = 0.1 sin(a;), (44) 

which corresponds to a growth rate u) — — 0.238J + 0.261. 

The computational domain has 200 equally spaced 
zones in the x-direction. We have repeated the run vary- 
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Figure 4. ID steady dynamo. Triangles represent cases with r\ = 
0.1 (and k = 1, 2, 5), squares cases with r\ = 0.05 [k = 1, 2), and 
diamonds with r] = 0.25 (fe = 1). Overplotted is the expected 
analytical solution with an exponential growth rate lui: the green 
dotted lines refers to cases with r] = 0.1, corresponding to Iu = 
0.385,0.567,0.59; the red dashed lines to cases with rj = 0.05, 
corresponding to Iu = 0.44, 0.77; and the solid blue line to r\ = 
0.25, corresponding to Iui = 0.236. 



ing 5 between 0.05 and 0.2, and the number of cells in the 
^-direction between 5 and 11, and found no appreciable dif- 
ference in the results. In Fig. 6 we show the evolution and 
compare the numerical results with the analytical expecta- 
tions. Both the drifting velocity of the wave and the expo- 
nential growth of its amplitude are properly recovered with 
an accuracy ~ 1%. Errors and convergence estimates are 
presented in Sect. 5.3 together with a comparison between 
the lst/2dn scheme and the IMEX scheme, and an estimate 
of the true growth rate and phase shift. 



5.2.3 Kinematic Couette flow 

This is a fully 2D dynamo test, on a non-Cartesian grid, cor- 
responding to a Couette flow. We use cylindrical coordinates 
R, z, with a domain extending in the radial direction in the 
range R = [0.1, 2], with 100 equally spaced zones, and along 
the 2-direction in the range z = [—1,1], with 100 equally 
spaced zones. This corresponds to an aspect ratio ~ 1. 

Differential rotation and density (and pressure) strati- 
fication are imposed such that the system is in equilibrium, 
by assuming a polytropic relation p oc p 4 ^ 3 and requiring 
that the Bernoulli integral 

In(h) + iln(l-w%)-^(r2-r2 ) 2 (45) 

is constant in the domain, where h — l + 4p/p is the specific 
enthalpy, = Q is the angular velocity, Qo is the angular 
velocity in R = 0, and the parameter A is an indicator of 
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Figure 5. ID steady dynamo, transition to the fastest growing 
mode. The left panel shows the B y component of the magnetic 
field, normalized to its maximum value. The transition from the 
k = 9 mode to the k = 5 mode happens around t = 12. The right 
panel shows the evolution of the maximum value of B v . The solid 
blue line is the expected growth for a k = 9 mode, the solid red 
line is the expected growth for the fastest k = 5 mode. 



the amount of differential rotation, since 

«%/(l - «%) = QA 2 {Q - Q ). (46) 

This approach is equivalent to the one used in models of 
differentially rotating NSs (Komatsu, Eriguchi & Hachisu 
1989; Font, Stergioulas & Kokkotas 2000; Bucciantini & Del 
Zanna 2011). Here the following values are adopted: Qo = 
0.3, p(R = 0) = p(R = 0) = 10, A 2 = 5, which correspond 
to i/(i? = 2) = 0.16. 

The system is followed to a final time t = 460 cor- 
responding to a dynamo period. Initially a purely toroidal 
magnetic field is imposed in the computational domain 

= 0.001 if s/(r- l) 2 + (Z-0.5) 2 < 0.25, 
B* = -0.001 if \/{r-l)' 2 + (Z + 0.5) 2 < 0.25, (47) 

corresponding to two distinct loops, with a zero net mag- 
netic field in the domain. We want to remind here that the 
shape of the initial magnetic field does not matter, because 
the system always selects the fastest growing eigenmode of 
the dynamo equation. The dynamo coefficient is £ = 0.12, 
while the resistivity is set to 7} = 0.03. Periodic boundary 
conditions are assumed in the Z— direction, while a perfect 
conductor with E = is assumed at the R = 0.1 and R = 2 
boundaries. In Fig 7 we show the result of the evolution, 
focussing on the B® component. It is evident from the top- 
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Figure 6. ID thin shear dynamo. The upper panel shows the 
B v component of the magnetic field, normalized to its maximum 
value, as a function of time and space. The solid line represent 
the analytical phase shift <x 0.261t. The lower panel shows the 
evolution of the maximum value of B y . The red solid line is the 
analytical solution. 



Figure 7. 2D Couette flow dynamo. The upper-left panel shows 
the B$ component of the magnetic field, in R = 1, normalized to 
its maximum value, as a function of time. The upper left panel is 
a contour plot of B^, in the domain at t = 40 (solid cyan lines) 
and t = 400 (dashed red lines), after exactly one wave-cycle. The 
lower panel shows the evolution of the maximum value of B^. 



right panel that an eigenmode is selected, with a constant 
drifting speed and an exponential growth. 

5.2.4 Dynamical NS dynamo in full GRMHD 

As a final test we present here the growth of an a 2 dy- 
namo (no differential rotation or meridional flow) in a neu- 
tron star, in the full dynamical general relativistic regime 
under the XCFC approximation, as currently employed in 
the X-ECHO code. The initial conditions are derived using 
the XNS code (Bucciantini & Del Zanna 2011), to which 
the reader is referred for a description of parameters char- 
acterizing the initial equilibrium configuration. The NS has 
a central density p c = 1.28 x 1CP 3 (in geometrized units 
c = G = Mq — 1). The equation of state used for the initial 
settings corresponds to the polytropic relation p — 100p 2 , 
however the system is evolved using an ideal gas EoS as 
done for all NS models in Bucciantini & Del Zanna (2011). 
The initial magnetic field is purely toroidal and its distribu- 
tion follows the barotropic law: 

ar 2 sin 2 8B 4, = K m (ar 2 sin 2 6ph) m , (48) 

with magnetic parameters K m = 10~ 4 and m = 1. r is the 
radius and 6 is the polar angle. 

The 2D simulation is performed in spherical-like coordi- 
nates for the conformal flat metric assuming axisymmetry. 
The computational domain is 9 — [— n,ir], with reflecting 



boundary conditions on the axis, and r = [0, 10], with dipole- 
like reflecting conditions at the center and zeroth-order ex- 
trapolation at the outer radius. The dipole-like conditions 
are chosen to allow a dipolar magnetic and electric field 
where physical quantities like the radial and components 
of the magnetic field do not vanish in r — > 0. The resistivity 
is set r\ = 0.05, while the dynamo coefficient is £ = 0.1, uni- 
form in the computational domain. The problem is followed 
to a final time t — 200. The time evolution of the metric is 
computed once every 100 steps of the MHD part. We want 
to stress here that these values have been chosen for the 
sake of having a fast numerical run, and have no physical 
significance. The scope here is to perform a test of the code 
in its full dynamical regime, and not to address a particular 
physical problem in realistic conditions. 

In Fig 8 we show the result of the evolution. The up- 
per panel shows the magnetic configuration that is obtained 
at the end of the run. The lower panel shows the growth 
of the maximum value of the poloidal magnetic field during 
the evolution. Having adopted a uniform £ term, spurious 
dynamo action is also present in the atmosphere surround- 
ing the NS. As a consequence a magnetic field develops and 
grows in the atmosphere too, despite it being completely un- 
magnetized at the beginning. Again, this is just due to the 
unphysical choice of a uniform dynamo term. In the figure, 
for the sake of clarity, we plot the poloidal field just inside 
the star. It is interesting to note that the growth rapidly 
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NS Dynar 





200 



Figure 8. Dynamical NS dynamo. The upper-panel shows the 
poloidal magnetic fieldlines. Colors represent the value of the 
component of the magnetic field in units 10 15 G. The thick 
solid line represents the surface of the NS. The lower panel shows 
the evolution of the maximum value of poloidal magnetic field 



approaches an exponential behavior, but that the slope (the 
growth rate) changes at around t — 140. This corresponds to 
the transition to a shorter wavelength. In particular, at the 
beginning the mode that is excited and grows corresponds 
to a radial wavelength of the order of twice the stellar ra- 
dius (given the initial conditions on B^). At t = 140 an 
eigenmode with a shorter radial wavelength of the order of 
the stellar radius emerges, corresponding to a faster growing 
mode. 



5.3 Convergence and compartison between the 
lst-2nd and IMEX schemes 

In this section we discuss the convergence properties of our 
scheme and compare the different time-stepping approaches. 
Despite the simplicity of our lst-2nd scheme, we do expect 
that, depending on the problem, its accuracy might be closer 
to 1st order, as opposed to the IMEX which should preserve 
second order accuracy in all cases. We want to stress here 



that the IMEX scheme we have implemented only acts on 
the solution of the MHD equations, and it is not integrated 
with the metric solver. However in the choice of the XCFC 
approach, we already assumed that the metric should be 
slowly varying in time with respect to the plasma. In this 
case we do expect that the order of the scheme should not 
be affected by the metric solver. 

We want to point out here that in order to evaluate the 
convergence of an algorithm, one needs a test problem with 
a smooth flow structure at all times including the initial 
conditions, and well posed boundary conditions that pre- 
serve the overall accuracy. If a reference solution, either an 
analytical solution or a high resolution run (such that nu- 
merical resistivity is smaller than the physical resistivity) 
for comparison is not available, one can compute the order 
of convergence using relative errors at different resolutions. 
In Ideal Relativistic MHD such a test was designed by Del 
Zanna, Bucciantini & Londrillo (2003b), using a large am- 
plitude, circularly polarized Alfven wave. 

In Resistive Relativistic MHD no such test has been 
presented yet. The current sheet problem of Sect. 5.1.1, ad- 
mits an exact solution only in the limiting case of infinite 
pressure. For the value of the pressure that we have used, in 
line with the previous existing literature, dynamical effects 
are present, and they dominate the deviations with respect 
to the solution Eq. 40. However, being a one dimensional 
problem, it is possible to use a high resolution run as a ref- 
erence solution. In Fig. 9 we show the relative error on one 
component of the magnetic field: 



Lx[B v {t = t)\ 



1 

N 



E 



\BV(t = i)-B? et (t = i)\\ 
max.[Bv(t = i)] 



(49) 



where the sum is done over the value of the quantity in 
each of the N cells of the computational domain, and the 
reference solution at time t is obtained interpolating a high 
resolution run with 1600 grid points. It is interesting to note 
that in this test problem the IMEX and our lst-2nd time 
stepping algorithm have similar performances, even if the 
IMEX has a smaller absolute error. This is because the solu- 
tion is slowly evolving in time and the displacement current 
is small compared to the conduction current, which is due 
to spatial gradients of the magnetic field. For this reason the 
overall order of the algorithm is dominated by the order of 
the explicit part of the solver. Moreover the order seems to 
be somewhat in between 2 and 3, as expected from that fact 
that we use a third order CENO spatial reconstruction. 

For the relativistic dynamo closure, as opposed to the 
resistive case, the ID steady dynamos admit analytical so- 
lutions. The other multidimensional tests we have presented 
lack a correct analytical solution and, being multidimen- 
sional, it is computationally prohibitive to run a high res- 
olution reference case, however as we will show we can use 
relative errors at different resolutions. In Fig. 10 we show 
the relative error on one component of the magnetic field, 
defined as in Eq. 49 for the ID steady dynamo with k = 1, 
r\ = 0.1, £ = 0.25. These values have been selected in order 
for the mode with k = 1 to be the fastest growing mode. 
ID steady dynamos are rapidly evolving in time, and with 
small spatial gradients. The displacement current dominates 
over the conduction current. The IMEX scheme performs as 
previously, with an order which is again between 2 and 3 
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Figure 9. Convergence test for the resistive current sheet prob- 
lem. Relative L\ errors on B v at t = 10, with respect to a high 
resolution run with 1600 grid points, in logarithmic scale, as func- 
tion of the number of grid points. Diamonds are for the lst-2nd 
scheme, triangles for the IMEX scheme. Dashed and dotted lines 
reprensent scaling respectively as —2.5 and —2.7. 



Figure 10. Convergence test for the ID steady dynamo problem, 
with k = 1, n = 0.1, £ = 0.25. Relative L\ errors on B y at t = 20, 
with respect to the analytical solution, in logarithmic scale, as 
function of the number of grid points (per mode wavelength). 
Diamonds are for the lst-2nd scheme, triangles for the IMEX 
scheme. Dashed and dotted lines reprensent scaling respectively 
as —2.3 and — 1. 



for the same reason as before. The lst-2nd scheme instead 
reduces to 1st order as expected. 

We have also attempted to estimate convergence of the 
solution and the order of accuracy of our algorithm for the 
thin shear layer dynamo, Sect. 5.2.2. We do not have an 
analytic solution available to use as a reference solution, 
and it is computationally prohibitive to run a high resolu- 
tion case, so we proceed comparing errors in runs a different 
resolutions. Moreover we do have an analytic approximated 
solution and we have an expectation for the functional form 
of the eigenmode. In particular we expect the various quan- 
tities to evolve as e~ Iijt . So if we assume that the unknown 
true solution changes in time according to this exponen- 
tial growth, and that the error of our numerical solutions 
scales with the number of grid points N as N~ p , we can 
fit simultaneously for the growth rate, phase shift, and ac- 
curacy of the scheme. The result is shown in Fig. 11. This 
is equivalent to estimate the order of accuracy by compar- 
ing relative errors. We find that both for the IMEX scheme 
and for our lst-2nd scheme the best fit estimate for oj is 
lu = -1(0.2383 ± 0.0001) + (0.2626 ± 0.0005), to be com- 
pared with the expectation of the approximated analytical 
solution uj = — 0.238/ + 0.261. The best fit for the order gives 
p = 2 ± 0.1 for the IMEX and p = 1.5 ± 0.1 for the lst-2nd 
scheme. We want to remember here that we use extrapo- 
lation for the boundary conditions in the z— direction, and 
the accuracy of the solution also depends on the boundary 
conditions that are used. 

In all of our tests, we have verified that, at the resolu- 
tion at which they were performed, the relative errors of the 
lst-2nd scheme are of order a few 10 . Unless one requires 
higher accuracy, or if the problem involves slowly varying 
solutions with strong spatial gradients, or for discontinu- 
ous solutions, given its easy implementation, we deem this 
approach satisfactory. The IMEX scheme in general shows 
better convergence already at 25 points per wavelength of 
the eigenmode, and should be the algorithm of choice for 
more demanding cases. 
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Figure 11. Convergence fit for the thin shear layer dynamo prob- 
lem. Relative L\ errors on B y at t = 24, with respect to the 
best fit reference solution, in logarithmic scale, as function of the 
number of grid points (per mode wavelength) in the x— direction. 
Diamonds are for the lst-2nd scheme, triangles for the IMEX 
scheme. Dashed and dotted lines reprensent scaling respectively 
as —2 and —1.5. 



6 CONCLUSION 

In this paper we have presented a fully covariant mean-field 
q— dynamo closure of the resistive relativistic MHD equa- 
tions, and shown how it can be implemented within a code 
for numerical 3 + 1 General Relativistic MHD. In particu- 
lar we have upgraded the Eulerian Conservative High-Order 
code for GRMHD, in its version for dynamical spacetimcs 
(Bucciantini & Del Zanna 2011). The X-ECHO scheme em- 
ployes a fully constrained method for Einstein's equations 
based on the extended conformally flat condition (XCFC), 
but this particular choice poses no constraint on the appli- 
cability of our dynamo closure, which is unchanged in the 
Cowling approximation (a static GR metric) as well as in 
other formulations of the Einstein equations. We have shown 
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that our implementation of a dynamo effect is straight- 
forward for any numerical scheme that has already been 
extended to include resistivity, since the novel generalized 
Ohm's law is very similar to that for resistive GRMHD and 
does not pose any additional numerical difficulty. This is 
true in the simple case proposed here, where an isotropic re- 
sistivity and an isotropic a-dynamo term (thus simple pro- 
portionality between the electric and magnetic fields in the 
comoving frame) have been considered. Far more complex 
closures have been developed for non-relativistic MHD, how- 
ever it is not obvious if and how fully covariant equivalent 
formulations could be found at all. We stress again that de- 
termining realistic values for the parameters that are used 
in the closure, or even finding an appropriate closure, is usu- 
ally non trivial, and often requires the use of mesoscale in- 
formations, and extrapolation of flow properties to small 
unresolved scales, ft rests to be proved that a mean field 
approach can achieve full resolution of microscopic physics 
on the macro-scale. 

This is the first numerical implementation of a fully 
covariant dynamo closure for relativistic MHD in the gen- 
eral dynamical regime. Instead of starting with the simpler 
Minkowskian case, we have directly proposed and applied 
our model to full GRMHD, exploiting the 3+1 formalism, 
in both static and evolving spacetimes. We have adopted 
a fully constrained strategy, to retain the general philoso- 
phy of the ECHO and X-ECHO schemes. The charge den- 
sity is derived from Gauss theorem, and it is not evolved as 
an independent quantity like in previous formulations (no 
appreciable differences are observed in the numerical tests 
available in the literature). The solenoidal condition on the 
magnetic field is preserved to machine accuracy using the 
Upwind Constrained Transport method on a staggered grid. 
Stiff terms due to small resistivity in the evolution equa- 
tion for the (spatial) electric field have been treated both 
with a simple implicit time-stepping procedure, which al- 
lows us, contrary to previous works, to obtain the ideal case 
of a perfectly conducting plasma simply by setting the re- 
sistivity coefficient to zero, and not just as a limit for small 
resistivity, and with a more roboust IMEX scheme. We have 
shown that depending on the problem, and in particular on 
the relative importance of spatial versus temporal gradients, 
a simple lst-order implicit scheme can give satisfactory re- 
sults, while the IMEX scheme tends to give more reliable 
performances independently of them. 

We have presented a set of standard tests, easily im- 
plementable and, where possible, we have compared the nu- 
merical results with the analytical expectations, confirming 
the robustness of the implementation. The majority of the 
tests have been performed in the so-called kinematic regime. 
This choice has a physical motivation: mean-field dynamo 
action arises from small-scale motions, that are almost al- 
ways strongly subsonic. Their kinetic energy is negligible 
with respect to the internal or gravitational energy of the 
system. Dynamo action is supposed to be quenched once 
the strength of the magnetic field reaches equipartition with 
the kinetic energy of the turbulent motions. As such, one 
expects that dynamo amplified magnetic fields cannot reach 
values high enough to affect the global dynamics. However, 
to verify the stability and robustness of the implementation 
in a more demanding regime, which must be fully dynamic 
and closer to a physical application, we have also presented 



a full dynamical case applied to a NS in GRMHD with a 
time- dependent metric. 

Our results show that it is possible to implement within 
codes for numerical relativity, a closure that in principle can 
allow one to model effects arising from dynamics at small 
scales, that would be prohibitive to follow in global simu- 
lations. The general idea of sub-gridding modeling effects, 
has been widely developed in the context of classical fluid 
dynamics, but is still in its infancy in the field of numer- 
ical relativity. There are several outstanding problems of 
relativistic fluid dynamics, ranging from the origin of rela- 
tivistic engines, to their characterization to the dissipative 
evolution of their magnetic fields, that involve extensive dy- 
namical ranges in space and time. The development of a 
clever sub-gridding approach offers an interesting possibility 
to investigate and model physical processes that otherwise 
would require a resolution that would be computationally 
prohibitive. This might lead to a different approach and to 
a rapid advancement in the field. 

We plan to apply our code for dynamo in GRMHD to 
both accreting disks around BHs and in proto-NSs. These 
two different environments are fundamentally related to the 
more promising engines for GRBs, and might have impor- 
tant implications for the general modeling of core-collapse 
supernovae. Strong magnetic fields have also been invoked 
for Short GRBs, which are commonly considered to be a pos- 
sible electromagnetic counterpart of binary mergers. Indeed, 
much of the MHD modeling done until now has focused on 
the large-scale properties, but it has been shown that the ex- 
pected magnetic configurations can be highly unstable, and 
a turbulent cascade is expected. Understanding if and under 
what conditions a mean-field dynamo can operate, what is 
its efficiency, and the geometry and topology of the result- 
ing field, might help to put stronger constraints onto the 
environment within which they are supposed to operate. 
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APPENDIX A: DYNAMO SOLUTION FOR A 
RELATIVISTIC THIN SHEAR LAYER 

We present here an approximated analytical solution for a 
simple relativistic shear layer (see problems in Sect. 5.2.1, 
5.2.2). It is the relativistic generalization of a problem pre- 
sented by Arlt & Riidiger (1999). Let us consider a shear 
layer, with an extension in the ^-direction much smaller than 
in the ^-direction, where we assume that all quantities do 
not depend on y. Within this layer the velocity has the fol- 
lowing profile: v x — v z = 0, v y = Sz, where S is the shear 
parameter, and does not change in time. 

We look for solutions of the form f(z)e II -" t ~ kx) , where 
I = \J— 1. Symmetry tells us that the electric and magnetic 
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fields must have the following parities: 

B x , E x oc 



B V ,B Z ,E V ,E Z 
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aiZ 
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a,iz 



a I(ujt — kx) 



(Al) 



where the coefficients en = Bf, Ef, B\ 



complex numbers, 
differ for the various quantities. We will look for a solution 
up to a z 2 order. 

First, the divergence- free condition on the magnetic 
field implies Bf = IkBf/2. We start from Eq. 32 and from 
the induction equation for the magnetic field, which yield 



7 - 1/2 a( 7 1/ "B l ) = e'^dsiaEu + e kl ml3 l B m ). (A2) 

In the special case of a flat metric in Cartesian coordinates 
we have a = 7 = 1, /3 l = 0, then 

£ (25 (E z z 2 + E Z )- 2BE) + iB^kSz 2 
+ AB y r) + 2-BoS + 2/r;^w + 2E% = 
i(2Bl +z 2 (2ESS + IBSk)) 

- 2I(z 2 (IBgS + Bfrjk + t]E z oj - IE Z ) 
+ B%r 1 k + EZ(r)Q.-I)) = Q 
(-2£(SV - 1) {Biz 2 + Bf) - B% n {k 2 z 2 - 2) 

- 2{-IBlr)k + z 2 {-Ir ] E? ) kS + ir)E y oj + 2 V E Z S- 
E\S 2 z 2 + Ef) + E y {I V oj - S 2 z 2 + 1))) = (A5) 
2E\ - IB^uj = (A6) 
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(2IE y k - 2IB z uj + k{2IE\ + B%lj)z 2 ) = 0. 



(A8) 



If we impose these relations to be satisfied both at the 0(1) 
and 0(z 2 ) order we can get a set of 4 linear equations for 
the variables B§ , Bq 1; B§: 

- IBScj + + r ? cj) 4 (4a 3 B l )»?S 3 + B^r){-I + rjuj) 2 

(2S 2 + k 2 (l + I V oj)) + 2a 2 r l S 2 {B? ) {-2 - 2I V oj)+ 
k{-2IB y r)S + B z (-I + r)u))) + 2Ia(-I + jjw) 
{B v (l + 2/770; + r, 2 (4S 2 - a; 2 )) + r,S 2 {2B z S+ 
B y (-r,k 2 -Iuj + V oj 2 )))) =0 (A9) 
l/(-J + r)u) 2 (-Z 2 B z S - + V uj)(2B v 1 ri - B y V k 2 

+ B Z S - IB y uj + B y r)u> 2 ) + C(Bo + IBS V lv+ 
Ik{B z + B y V S + IBfau)) = (A10) 
B z u + (Jfc(£B» + r){B x + IB z k)))/{-I + rflj) = (All) 
B\lo + l/(2(-J + r 1 ujf)k{-2Ia i B! ) S 2 - 2{B y i)k 

+ IBSS)(-I + r,uj) 2 - 2^ 2 S(B y V kS + BZ(-I + V lu)) 

+ a(-I + r)u)(2S(2B y eta + B Z S) + B%k(-I + tjoj)))) = 0. 

(A12) 

Imposing that the determinant of the matrix representative 
of this system vanishes, provides us with a fifth-order equa- 
tion for uj as a function of k, r\, £, S which is our dispersion 
relation. 

For example, in the case £ = 0.2, r/ = 0.1, k = 1, 
and S — 0.9 we find a growing mode solution with u = 
0.261 — 0.2387. To this mode it corresponds an eigenvector 
with IBqI/IBqI = 2.13 and with a phase difference between 



these two components 8(f) = 0.66. The non-relativistic case, 
instead, where displacement currents and charge densities 
are neglected, and where an exact analytical solution can be 
found, gives u = 0.254 - 0.240J. 

In the case S = it is possible to find an exact analytical 
solution for the dispersion relation and the eigenmodes 



- r,k 2 

which gives: 



Ioj + t)oj 2 )/(-I + -nuj) = 0, 



u = /[l ± v /l + 477fc(£-r?fc)]/(2 J y). 



(A13) 



(A14) 



Exponentially growing modes are possible only for k < £/r). 
The quantity £,/"qk is the dynamo number. The fastest grow- 
ing mode has k = £/{2rj) and a growth rate oj max = 
(a/I + £ 2 — l)/(2r?). The eigenmode corresponding to the 
growing solution is 



B y = IB Z ,B X = 0,E X = 



E y = 



E z 
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This solution can be compared with the non-relativistic case, 
where the displacement current is neglected 



uj = Ik(a — r/k) 
B y = IB Z ,B X = 0. 
E y = £B y + IB z r,k 
E z = £B Z - IB y r)k. 



(A16) 
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